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W). Abstract 
c/3 ' 

^ ' The time evolution of the strength of the Earth's virtual axial dipole moment 

* ^ I (VADM) is analyzed by relating it to the Fokker-Planck equation, which describes a 

random walk with VADM-dependent drift and diffusion coefficients. We demonstrate 
[ first that our method is able to retrieve the correct shape of the drift and diffusion 

, coefficients from a time series generated by a test model. Analysis of the Sint-2000 data 

shows that the geomagnetic dipole mode has a linear growth time of 20^7'^ kyr, and 
that the nonlinear quenching of the growth rate follows a quadratic function of the type 
[ [1 — {x / xq)'^] . On theoretical grounds, the diffusive motion of the VADM is expected to 

' be driven by multiplicative noise, and the corresponding diffusion coefficient to scale 

I quadratically with dipole strength. However, analysis of the Sint-2000 VADM data 

' reveals a diffusion which depends only very weakly on the dipole strength. This may 

I I indicate that the magnetic field quenches the amplitude of the turbulent velocity in 

' the Earth's outer core. 

I> 

, Keywords: Geodynamo, Reversals, Secular variation, Sint-2000 record, Turbulent con- 

vection. Stochastic processes. 
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^ ! 1 Introduction 



The strength of the geomagnetic dipole moment shows a considerable time variability, 
about 25% r.m.s. of the mean, over the course of thousands of years. Occasionally, the 
variability is so large that the sign of the dipole moment changes. These reversals happen 
roughly once per (2 — 3) x 10^ yr (Merrill et al., 1996). The geomagnetic field is the 
result of inductive processes in the Earth's liquid metallic outer core. Helical convection 
amplifies the magnetic field and balances resistive decay. Several groups have confirmed 
this idea with the help of numerical simulations (Glatzmaier and Roberts, 1995; Kuang 
and Bloxham, 1997; Christensen et al., 1999). A suitable measure of the geomagnetic 
dipole is the Virtual Axial Dipole Moment (VADM), of which several records have been 
published, e.g. by Guyodo and Valet (1999) and Valet et al. (2005). Since the dipole 
moment is the result of many processes taking place in the convecting metallic outer core 
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Figure 1: The Sint-2000 VADM data of Valet et al. (2005), a time series of 2000 unsigned VADM 
values spaced in time by 1000 year, covering the past 2 Myr history. We have inserted a sign flip 
at the times of known reversals, see text. 



that interact with each other in a complicated way, it makes sense to try to describe the 
time evolution of the VADM over long time scales as a stochastic process. 

Before entering into details we recall that statistical modelling of the geomagnetic field 
has a long history. Constable and Parker (1988) were the first to give a complete char- 
acterization of the statistical properties of the geomagnetic field in terms of its spherical 
harmonic expansion coefficients. The distribution of the axial dipole was found to be sym- 
metric and bi-modal, consisting of two gaussians shifted to the peak position of the two 
polarity states. They also showed that the expansion coefficients of the non-dipole field 
may, after appropriate scaling, be regarded as statistically independent samples of one 
single normal distribution with zero mean. This GGP (giant gaussian process) approach 
as it is now generally referred to, permitted computation of the average of any field-related 
quantity. Hulot and Le Mouel (1994) have extended the GGP approach by considering the 
evolution of the statistical properties with time, and Bouligand et al. (2005) have tested 
the GGP modelling technique on hydromagnetic geodynamo simulations. 

Returning to the time evolution of geomagnetic dipole as a stochastic process, consider 
a stochastic equation of the type 

X = v{x) + F{x)L{t) . (1) 

The function v{x) has the dimension of a velocity and represents the effective growth 
rate of x, sometimes called the drift velocity. The fluctuations are embodied in the term 
F{x)L{t) and they induce an additional diffusive motion of x0 Here L{t) is a stationary 
random function with zero mean and a short correlation time Tc'- 

{Lit)) = 0, (L(t)L(t - r)) = l2^,3 Te 5(t) . (2) 

A short correlation time means that the duration Tc of the memory of L{t) is much shorter 
than all other time scales in the process. Under these circumstances the autocorrelation 
function of L{t) behaves as a (5-function of time. The probability distribution p{x, t) of x{t) 
determined by Eq. ([1]) obeys the Fokker-Planck equation (Van Kampen, 1992; Gardiner, 
1990):| 

Here t is time, and v is again the effective growth rate of x. The diffusion coefficient is 
equal to 

POO 

D ^2F^ {L{t)L{t -T))dT ^ F^L^^ ^ , . (4) 
Jo 



^The noise is called additive if F is constant, otherwise it is referred to as multiplicative noise. 
^Provided vtc ^ x; this particular form of Eq. Q requires in addition that dD/dx <^ v. 
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Figure 2: The amplitude distribution of the unsigned Sint-2000 data. In principle the distribution 
is symmetric with respect to VADM = 0, and has a characteristic double-hump structure with 
small but nonzero probability at VADM — 0. 



The Fokker-Planck equation is a simple and versatile tool for modelling the dynamics of 
a stochastic process. That is to say, the statistical properties of a wide variety of different 
stochastic processes can be described by the Fokker-Planck equation ([3]). Hoyng et al. 
(2002) have shown that for theoretically plausible functions v{x) and D{x) the amplitude 
distribution of the Sint-800 data (Guyodo and Valet, 1999) is very well predicted by 
Eq. (H. 

The purpose of this paper is twofold. We investigate whether the Sint-2000 VADM 
time series (Valet et al., 2005) can indeed be described by a Fokker-Planck equation 
([3]). Secondly, we derive the dependence of the effective growth rate v and the diffusion 
coefficient D on the magnitude x of the VADM without making any prior assumption on 
the functional form. In doing so we are able to measure the linear growth rate of the 
dipole mode and its nonlinear quenching from the data. Likewise, the diffusion coefficient 
D{x) provides information on the convective flows in the outer core. This marks the 
difference between our approach and that of the GGP: we do not stop at giving a statistical 
desciption of the multipole coefficients of the geomagnetic field, but we extract information 
immediately related to the physics of the geomagnetic dipole. 

After a brief discussion of the Sint-2000 data in Section 2, we develop in Section 3 a 
technique for extracting the functions v{x) and D{x) from a time series. Next, in Section 
4, we validate the method with the help of an artificial VADM time series generated by a 
simple model to see how well we can retrieve the v{x) and D{x) that were used to generate 
the series. In Section 5 we apply the method the Sint-2000 VADM data (Valet et al., 2005) 
and we discuss the implications of our findings for the geodynamo. A summary and our 
conclusions appear in Section 6. 

2 Sint-2000 data 

The Sint-2000 data comprises a time series of 2000 unsigned VADM values spaced by 1000 
year, covering the past 2 Myr history of the geomagnetic dipole. The positions of the 
reversals are indicated in Figure 2 of Valet et al. (2005), where they show up as local 
minima in the VADM record. To obtain a VADM time series with sign we have inserted 
a sign change between those locations. The result is shown in Figure [H In doing so we 
may miss some of the fine structure in the reversal time profile. However, in view of the 
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considerable intrinsic uncertainties in the data [see Figure 2 of Valet et al. (2005)], similar 
and larger ambiguities apply to the whole VADM time series. The amplitude distribution 
of the unsigned VADM data is shown in Figure [2j 



3 Method 

We start with the discretized version of the Fokker-Planck equation. Discretization of 
space and time in Eq. [3] leads to 



Pi{t + At) - pi{t) _ Vi+ipi+i{t) - Vi-ipi-i{t) 

(5) 



At 2Ax 

Di+ipi+i{t) + A-i/Oi-i(*) - 2DiPi{t) 



2(Ax)2 

We may rewrite this equation in matrix form as 

Pi (t + At ) = ^ {5ij + AtMij ) pj (t) , (6) 
j 

where is M a tridiagonal matrix with elements 



2Ax 2(Ax)2 



1,1 



(Ax)2 ' 



- "2A^ + 2(A^- 

By solving for Vi and Di we obtain 

Vi = (M,+i,i - Mi_i,i) Ax ; 

Di = (M,+i,, + Mi_i,i) (Ax)2 . (8) 

These expressions will be used to infer Vi and Di once the matrix M has been determined. 

Each column of M adds up to zero, Mj_i^j + Mj.j + Mj+i^j = 0, so that there are two 
free parameters per spatial interval i, exactly as many as the parameters Vi and Di in 
the Fokker-Planck equation. An important consequence of the zero column sum is that 
Eq. ([UD is norm-conserving, 

Y,p,{t + At)=Y,p,{t) . (9) 



In the limit At — Eq. ^ becomes 

dpi 
dt 



^Mi^pj. (10) 



The object of this paper is to extract the effective position-dependent (= VADM-dependent) 
velocity and diffusion coefficient from a time series, in this case of the strength of the 
Earth's magnetic dipole moment. To this end we construct from the data the matrix T 
whose elements Tij contain the transition probabilities for a system in position j at some 
time t to move to position i at a later time t + r. We begin by counting the number of 
times Nij that the system is located in position j at some time t and in position i at time 
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Figure 3: The left panel is the transition matrix T as obtained from simulation data of the HD 
model with r = 0.007. Right panel: the approximate transition matrix T = exp(TAf), where M is 
a tridiagonal matrix, see text. The upper left corner corresponds to (—2,-2), the lower right to 
(2, 2). The bin size is 0.16 x 0.16. The matrix elements obey < Ty < 1 and J2t Tt, = 1. 



t + T. The required matrix elements are then equal to Tij = aj ■ Nij and have a statistical 
error aij = Uj ■ Nij^/"^ . The normalization coefficients Oj are fixed by the requirement 
that the columns of Tij should add up to unity, Tij = 1. The time lag r, finally, must 
be chosen comparable to, or larger than the correlation time of the randomly fluctuating 
part of the system, but small in comparison to the time scale on which the data changes 
systematically. 

Our assumption is that the process is Markovian and therefore can be described by 
Eq. (fTO|) . from which it follows that pi{t + t) = J2j^w{'''M)ijPj{'^)- The theoretical 
transition matrix T is therefore 

T = exp(rM) , (11) 

and our goal is now to find a tridiagonal matrix M such that T closely resembles T. The 
matrix M has approximately 3n degrees of freedom (ignoring boundary effects) , of which 
n can be eliminated by norm conservation (columns add up to zero). To find the remaining 
2n degrees of freedom we minimize the function 




exp(TM)i 



CTij 



(12) 



We could follow an alternative approach, by using the stationary distribution pi = pi{oo) 
which we may find by binning the data as in Figure [21 Since the stationary distribution 
should obey Eq. ([TO]) , we have J2j ^ijPj = 0- This relation can be used to eliminate 
another n degrees of freedom in M, after which we find the remaining n by minimizing 
the function (fT2]) . But we opted for fitting 2n degrees of freedom and to use J2j ^ijPj — 0) 
or equivalently X]j '^ijPj = , as a consistency check on our computations. 
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Figure 4: The diffusion coefficients Di as a function of x, obtained by fitting the simulation data 
of Hoyng and Duistermaat (2004) to the Fokker-Planck equation. The drawn line is given by (fT4|) 
with Do = 0.4 and (7-^) = 0.27. The inset shows the effective growth rate Vi, compared to the 
theoretical value x{l — x'^) (drawn line). The error bars indicate 80% confidence intervals. 

4 Validation with the HD model 

First, we test the approach outlined in the previous section on data generated with the 
model of Hoyng and Duistermaat (2004) i This is a time series x{t) of VADMs measured 
in units of the equilibrium value, so x = 1 corresponds to the nonlinear equilibrium value 
of the VADM. The series comprises 5 x 10^ data points with a time spacing of 0.001, 
and by construction this is also the correlation time0 Time is measured in units of the 
linear growth time of the dipole mode so that the series is about 50 Myr long in real time. 
We discretize the strength x(t) of the VADM into 25 bins of width 0.16 in dimensionless 
units, and we construct a histogram of all sets {x(t),x{t + r)} employing a time lag of 
r = 0.007. We then exploit the fact that there is no sign preference, that is, for a given 
realisation x{t) the series —x{t) is an equally likely realisation. Accordingly, we add to the 
histogram all sets {—x{t), — x(t + r)}. We follow the procedure of the previous section, and 
the resulting effective transition matrix T is plotted in Figure [3l left panel. Note that the 
blue matrix elements near the centre have a relatively large value but not a more accurate 
one: matrix elements near the centre of the figure are determined by small x{t) associated 
with reversals and these are rare. The most accurate elements correspond therefore to 
the equilibrium value x = ±1, and are located in the wings near (1,1) and (—1,-1) in 
Figure [3] (left panel) . 

We then perform the fitting procedure outlined above to obtain the tridiagonal transi- 
tion matrix M, and in the right panel of Figure [3] we have plotted T = exp(rM). There 
is a clear similarity between the two matrices. Figure |4] shows the resulting values for the 

^Henceforth referred to as 'HD' or 'HD model'. 

*The time resolution of this series is a factor 10 higher than that of the series used in Fig. 2 of HD, 
but the other parameters are the same {a = 2, c — 5, and D = 0.4). 
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diffusion coefficient D and the effective growth rate v (inset), computed with the help of 
([8]). The error bars are 80% confidence intervals computed with the bootstrap method 
(Newman and Barkema, 1999). This captures the statistical errors, but not the systematic 
errors. 

4.1 Comparison with theory 

To place these results in perspective, we compute the Fokker-Planck equation for the 
probability density of x by integrating Eqs. (5) and (6) of HD over the overtone amplitude 
r, to find 

| = _^,(:_,>+i|,o„,.^ + (.^,U)p. (13) 

Hence, we recover Eq. ([3]) with 

v = x{l-x^); D = Do{x^ + {r^)\^) . (14) 

Here Dq a constant equal to 0.4 for the HD dataset used here, and {r'^)\x is the mean square 
overtone amplitude for given x. The result is the two drawn lines in Figure [H Since {r'^)\x 
is only a weak function of x, we did not bother to measure it from the simulation data. 
Instead, we replaced it by the average of over all x, measured to be (r^) = 0.27. The 
V and D recovered from the data compare rather well with their theoretical values (I14p . 
The agreement for D could be further improved by allowing for the fact that {r'^)\x is 
smaller than 0.27 near x = and larger than 0.27 for x > 1. However, we cannot expect 
agreement to within the statistical errors because of approximations made in deriving the 
Fokker-Planck equation (1131) . As a result there are small systematic differences between 
the statistical properties of x{t) predicted by ^T3\\ and (fH|) and those of the numerically 
generated x{t). These differences are visible because we use many data points (5 x 10^). 

These results demonstrate that our analysis is capable to extract the information on 
the effective VADM growth rate v and the type of noise that was used to generate the 
time series. The scaling D oc is a consequence of the multiplicative noise that the HD 
model employs [that is, a noise term of the type x = • • • + N{t)x]. But that is really a 
detail here. The main issue is that we have succesfully validated our retrieval method, as 
we have shown that our analysis is able to get out what has been put into the model. To 
avoid misunderstanding we note that this agreement does not say anything on whether 
the HD model describes the physics of the geomagnetic dipole correctly or not, or better 
than other reversal models do. It only tells us that our retrieval method appears to work 
satisfactorily. 

5 Application to the Sint-2000 data 

We then repeat the same procedure on the Sint-2000 data. Figured) The fitting procedure 
was performed with a time lag of r = 4 kyr, see Figure [5l This choice is motivated as 
follows. The autocorrelation time of VADM data is a few hunderd yr (the time scale for 
rapid random changes in the geomagnetic dipole), but the sampling of the Sint-2000 VADM 
data increases that to 1 kyr. The time scale for systematic changes may be identified with 
the linear growth time of the dipole mode (of the order of 10 kyr). The resulting effective 
growth rate v and diffusion coefficient D are shown in Figure [6l The error bars are again 
80% confidence intervals computed with the bootstrap method (Newman and Barkema, 
1999). In reality, the errors will be larger as we did not allow for the considerable intrinsic 
errors in the VADM data (Valet et al., 2005). 
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Figure 5: Left panel: the transition matrix T obtained from the Sint-2000 data, with t — 4 kyr. 
Right panel: the approximate transition matrix T — exp(MT) obtained from a tridiagonal matrix 
M. The upper left corner corresponds to (—10,-10) • 10^^ Am^, and the lower right corner to 
(10, 10) • 10^^ Am^. We use square bins of linear size 1 • 10^^ Am^. 

The X dependence of the effective growth rate is approximately as expected. The best 
fit of the function Ax[l — (x/xq)^] to the 'data points' Vi yields 1/A = 20^*^7^ ky, and 
xq = (5.4 lb 0.5) X 10^^ Am^. For small x we have v oc x which corresponds to linear 
growth of the dipole mode when it is small, and the —x^ term is the nonlinear quenching. 
The surprise is in the x-dependence of the diffusion coefficient D which we discuss below. 

5.1 Implications for the geodynamo 

The analysis of the Sint-2000 data confirms that the geomagnetic dipole mode is unstable 
with a linear growth time 1/A ~ 20^^^ ky. The nonlinear quenching follows approximately 
a quadratic quenching function [1 — (x/xq)^]. The nonlinear equilibrium is attained at a 
VADMofxo = 5.4-10^^ Am^. These results are more or less as expected. To our knowledge 
this is the first time that the linear growth rate A and the shape of the quenching function 
of the geomagnetic dipole have been measured from pertinent data. 

In order to judge our results on the diffusion coefficient we derive the theoretical x- 
dependence of D. To this end we consider the induction equation of MHD: dB/dt = 
V X (vq X B) + V X (6v X B) + 7]V'^B. The fluctuating velocity 5v represents the convective 
turbulence in the metallic outer core superposed on a steady flow vq. If we expand the 
magnetic field in the induction equation in multipoles, the equation for the dipole becomes 
X = • • • + const • 5v{t)x. Only the contribution of the fluctuating term acting on the dipole 
is written down explicitly. Comparing with Eq. ([T|) which produces a diffusion coefficient 
D = 2F'^ J^{L{t)L[t — T))dT, we now obtain D oc ((5f )^ ^ Tc x^ cx /3x^, where Tc is 
the correlation time of 6v{t), and /3 ~ ((^f )r.m.s.'^c the turbulent diffusion coefficient that 
occurs in the dynamo equation (Moffatt, 1978). Detailed considerations lead to (Hoyng, 
2008) 

P x2 

^ ^ E^IV ^ const , (15) 

where R is the radius of the outer core and the number of convective cells in the 
core. There is a small, approximately constant contribution to D due to feedback of the 
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Figure 6: Diffusion coefficients Di and velocity Vi (inset) as a function of magnetic dipole strength, 
obtained from fitting the Fokker-Planck equation to the Sint-2000 data. The drawn line in the 
inset is the best fit of A2;[l — (x/xq)^] to the 'data points' Vi, see text for details. 

overtones on the dipole amplitude. This term also occurred in the HD model, cf. Eq. (jl4p . 
It is important because it is related to the occurrence of reversals but it plays no role in 
the following discussion. 

We expect therefore that D cc x^, and the explanation is simple. The form of the 
induction equation makes that a given 5v generates a change in B proportional to the 
magnitude of B. For given 5v the diffusive motion of B is therefore larger if B is large, 
and this translates to the dipole component as well. However, these considerations are 
not borne out by our numerical results in Figure [6l 

The increase of the diffusion coefficient for |x| in Figure [6] is probably an artifact 
of the restricted length of the data, in combination with the fact that there are only five 
reversals and one aborted reversal in the last 2 Myr. VADMs smaller than 2 x 10^^ Am^ 
are absent in the Sint-2000 data except during the very brief reversal periods. Since the 
data cannot resolve the fine structure of the VADM during a reversal one might wonder 
what the effect would be of a few rapid sign changes near a reversal, but that would only 
serve to make -D(O) larger. 

For VADM > 2 x 10^^ Am^ the diffusion coefficient is constant. In fact, one might say 
that the data are consistent with a constant D at all VADM. We have tested the possibility 
that this result might somehow be caused by the limited time resolution of the Sint-2000 
data. To this end we have generated from the HD data a Sint-2000-like series by taking 
a running average and then a subset of 2000 data points separated by 1000 yr. This can 
only be done approximately as we cannot convert the dimensionless time of the HD data 
(in units of the dipole growth time) into real time. It is conceivable that this new time 
series would have a large D{0) and a constant D at larger x, but the resulting v and D 
did not differ materially from those in Figure [H 

There seems to be no signature of multiplicative noise in the Sint-2000 data, and we 
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believe that this is a sohd result. Instead, the data indicate that the noise is quasi- 
additive, because the diffusion coefficient D ~ F'^L'^ ^ ^ t^ in Eq. ^ would be constant 
if F in Eq. ([T]) is constant. We have no explanation for this, but there is the intriguing 
possibility that it is due to a nonlinear quenching of the fluid velocity fluctuations 6v{t). 
If /3 oc {{dv)'^)Tc would scale as cx l/B"^ oc then D would be effectively independent 

of VADM, cf. (fT5]) . The best way to test these ideas would be to use a longer dataset, 
and an obvious possibility would be to use VADM data from hydromagnetic geodynamo 
simulations for this purpose. 

6 Summary and conclusions 

We have presented and validated a technique for extracting the effective growth rate and 
diffusion coefficient of a time series of a stochastic process. An attractive feature of the 
method is that it does not assume any a priori mathematical form for these quantities. 
Application of the method to the Sint-2000 VADM time series has shown that it is possible 
to measure the linear growth rate of the geomagnetic dipole and the shape of the nonlinear 
quenching of this growth rate. The dependence of the diffusion coefficient on the VADM 
suggests that the amplitude of the convective flows in the outer core is suppressed with 
increasing dipole strength. The main limitation in extracting more useful information on 
the geodynamo is the length of the Sint-2000 series. In future research, we will apply this 
analysis technique to time series obtained from hydromagnetic geodynamo simulations. If 
no nonlinear quenching is observed in these simulations, the simulation model produces 
time series which are qualitatively different from the Sint-2000 VADM time series. On the 
other hand, if these simulations show similar nonlinear quenching, then the cause of it can 
be investigated within the model. 
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